!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the program is meant to design the wall of a supersonic nozzle, 
! which provides a maximum thrust at exit with constraints of 
! length and ambient pressure
! ref to: 1. chapter 16 of 'Gas Dynamics' by Maurice J. Zucrow & 
!            Joe D. Hoffmann
!         2. G. V. R. Rao. Exhaust nozzle contour for optimum thrust
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-07-24
! revised  : 2015-07-29 add thetaB loop and re-code it
!            2015-08-05 specify length and ambient pressure (more real)
!            2015-08-09 re-calcuate thrust on the transonic line use p0e before end
!            2015-08-10 add quadratic curve fitting in region R calculation to fining
!                       mass flow rate determination
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! called subroutine:
!   VariableDef (module)
!   Thermo
!   Transonic
!   InnerP
!   SymmetryP
!   InverseW
!   Vtheta
!   Thrust
!   Interp
!    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   i,j    counter
!   k~k5   counter
!   kd     acounter for point D along prescribed c-line
!   flagd,flagd0 if current error product greater than previous one, this theta and 
!          previous theta loops, 0-no
!   k20    point D number on previous theta loop
!   kinterp point number where interpolation included along the previous r-c-line
!   ffn    array to store every c-line point amount
!   xbarma array to store the axial location of Mach one line
!   mdotDE0 mass flow rate along DE, from D to E-1 
!   atemp~ytemp temp variables of a M alpha y
!   interc interpolate coefficient
!   thetai angle step in inverse wall region to determine the point along the TB
!   stepd  step length in y used to distribute the exit line DE
!   thetad streamline angle respect to x-axis at point D
!   alphad Mach angle at point D
!   C2~C3  lagrange mutiplier coefficients
!   p0e	  p0 calculated from pe
!   mul*   error product, e- E loop, d-D loop, t-t loop, 0-previous loop
!   mdot1  mass flow rate through any region ending c-line
!   arou~ap average variables
!   mdotBD mass flow rate through line BD
!   mdotDE mass flow rate through line DE
!   dmdot  mass flow rate used in region R calcualtion
!   dmdot0 mass flow rate before previous point in region R calcualtion
!   mdotDF mass flow rate through DF
!   x~M    arraies of point properties
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning :
!   thetaB must not too large, for examplexed < 5.0 and Med > 2.8
!   p0 must not set too small because of balance of error ans time
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program MAXT_IRRO()
   use VariableDef
   implicit none
!   temporary variables
   integer::i,j,k,k1,k2,k3,k4,k5,kd,flagd0,flagd,k20,kinterp
   integer,allocatable::ffn(:)
   real*8,allocatable::xbarma(:),mdotDE0(:)
   real*8::atemp,Mtemp,alptemp,ytemp,interc
   real*8::thetai,stepd,thetad,alphad,C2,C3,p0e
   real*8::mule0,mule,muld0,muld,mult0,mult
   real*8::mdot1,arou,aV,aVu,ay,atheta,aalpha,ap,mdotBD,mdotDE,dmdot,dmdot0,mdotDF
!   output variables
   real*8,allocatable::x(:,:),y(:,:),Vu(:,:),Vv(:,:),V(:,:),p(:,:),T(:,:)
   real*8,allocatable::rou(:,:),M(:,:)
   
!   **************************************** program start *************************!
!**************************** initiate variables **************************
!****************    Initiate variable     ***

! initiate the variables
   el=0.0001
   ev=0.1
   pi=3.141592653
   delta=1
 
   Tt=3000
   Pt=7000000
   g=1.2
   Rg=320.0
   p0=50000.0
   
   yt=1.0 !0.1 !0.025 !1.0
   rtu=2.0 !0.2 !0.050 !2.0
   rtd=0.5 !0.05
   thetai=0.5*pi/180
   stepd=yt*0.025
   xed=8.0
   
   wallflag=1
   flagy40=0
   
   Ny=11
   Ni=299
   Nj=299
! allocate memory for array
   allocate(xbarma(Ny))
   allocate(ffn(Ni))
   allocate(mdotDE0(Ni))
   
   allocate(x(Ni,Nj))
   allocate(y(Ni,Nj))
   allocate(Vu(Ni,Nj))
   allocate(Vv(Ni,Nj))
   allocate(V(Ni,Nj))
   allocate(p(Ni,Nj))
   allocate(T(Ni,Nj))
   allocate(rou(Ni,Nj))
   allocate(M(Ni,Nj))

!***********************  initial line ******
   call Transonic(xbarma,x(1:Ny,1),y(1:Ny,1),Vu(1:Ny,1),Vv(1:Ny,1),V(1:Ny,1)&
         &,M(1:Ny,1),T(1:Ny,1),p(1:Ny,1),rou(1:Ny,1))
   ffn(1)=1

!**********************   initial region *********
   do k1=2,Ny
      do k2=2,2*k1-2
         x1=x(k1,k2-1)
         y1=y(k1,k2-1)
         Vu1=Vu(k1,k2-1)
         Vv1=Vv(k1,k2-1)
         x2=x(k1-1,k2-1)
         y2=y(k1-1,k2-1)
         Vu2=Vu(k1-1,k2-1)
         Vv2=Vv(k1-1,k2-1)
         
         call InnerP()
         if (flagy40==1) then
            flagy40=0
            exit
         end if
         
         x(k1,k2)=x4
         y(k1,k2)=y4
         Vu(k1,k2)=Vu4
         Vv(k1,k2)=Vv4
         V(k1,k2)=V4
         M(k1,k2)=M4
         T(k1,k2)=T4
         p(k1,k2)=p4
         rou(k1,k2)=rou4
      end do
      x1=x(k1,k2-1)
      y1=y(k1,k2-1)
      Vu1=Vu(k1,k2-1)
      Vv1=Vv(k1,k2-1)
      
      call SymmetryP()
      
      x(k1,k2)=x4
      y(k1,k2)=y4
      Vu(k1,k2)=Vu4
      Vv(k1,k2)=Vv4
      V(k1,k2)=V4
      M(k1,k2)=M4
      T(k1,k2)=T4
      p(k1,k2)=p4
      rou(k1,k2)=rou4
      ffn(k1)=k2
   end do
!   check mass flow rate balance, right
   mdot1=0.0
   do i=2,2*Ny-1
      arou=(rou(k1-1,i-1)+rou(k1-1,i))/2
      aV=(V(k1-1,i-1)+V(k1-1,i))/2
      atheta=(atan(Vv(k1-1,i-1)/Vu(k1-1,i-1))+atan(Vv(k1-1,i)/Vu(k1-1,i)))/2
      aalpha=(asin(1/M(k1-1,i-1))+asin(1/M(k1-1,i)))/2
      ay=(y(k1-1,i-1)+y(k1-1,i))/2
      mdot1=mdot1+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)&
            &*(y(k1-1,i-1)-y(k1-1,i))
   end do

!****************************** inverse wall region and point D location  ************

   thetab=0.0   ! initiate thetaa
   mult0=1e6   ! initiate mult0 with a huge number
   flagd0=0   ! initiate flagd with 0 in previous theta loop
   k20=ffn(Ny)   ! initiate k20
   k=2*Ny-1   ! initiate the first inverse wall c-line point number
   do   ! incresing thetab along arc TB, inverse wall, expansion wall
      thetab=thetab+thetai	!thetaB increased by step of thetai
      ! calculate whole inverse wall right c-line
      ! there the next c+ characteristics line must not intersect the wall
         !before the prescribed point
      ! the l3 is the number of which the front c-line point start intercept
         !current r-c-line by l-c-line 
      do k3=2,k   ! k3 loop begin
         x(k1,1)=rtd*sin(thetab)
         y(k1,1)=yt+rtd*(1-cos(thetab))
         xb=x(k1,1)   ! store point B location
         yb=y(k1,1)
         
         x3=x(k1-1,1)
         y3=y(k1-1,1)
         Vu3=Vu(k1-1,1)
         Vv3=Vv(k1-1,1)
         x4=x(k1,1)
         y4=y(k1,1)
         
         call Thermo(V(k1-1,k3),p(k1,1),rou(k1,1),T(k1,1),atemp,Mtemp)
         alptemp=atan(Vv(k1-1,k3)/Vu(k1-1,k3))+asin(1/Mtemp)
         ytemp=y(k1-1,k3)+(x(k1,1)-x(k1-1,k3))*tan(alptemp)
         
         if (ytemp < y(k1,1)) then
            x1=x(k1-1,k3)
            y1=y(k1-1,k3)
            Vu1=Vu(k1-1,k3)
            Vv1=Vv(k1-1,k3)
            
            call InverseW()
            
            Vu(k1,1)=Vu4
            Vv(k1,1)=Vv4
            V(k1,1)=V4
            M(k1,1)=M4
            T(k1,1)=T4
            p(k1,1)=p4
            rou(k1,1)=rou4
            exit
         end if
      end do   ! k3 loop begin
      
      ! successive point on the r-c-line
      do k2=2,k-k3+2
         x1=x(k1,k2-1)
         y1=y(k1,k2-1)
         Vu1=Vu(k1,k2-1)
         Vv1=Vv(k1,k2-1)
         x2=x(k1-1,k2+k3-2)
         y2=y(k1-1,k2+k3-2)
         Vu2=Vu(k1-1,k2+k3-2)
         Vv2=Vv(k1-1,k2+k3-2)
         
         call InnerP()
         if (flagy40==1) then
            flagy40=0
            exit
         end if
         
         x(k1,k2)=x4
         y(k1,k2)=y4
         Vu(k1,k2)=Vu4
         Vv(k1,k2)=Vv4
         V(k1,k2)=V4
         M(k1,k2)=M4
         T(k1,k2)=T4
         p(k1,k2)=p4
         rou(k1,k2)=rou4
      end do
!      axial point
      x1=x(k1,k2-1)
      y1=y(k1,k2-1)
      Vu1=Vu(k1,k2-1)
      Vv1=Vv(k1,k2-1)
      
      call SymmetryP()
      
      x(k1,k2)=x4
      y(k1,k2)=y4
      Vu(k1,k2)=Vu4
      Vv(k1,k2)=Vv4
      V(k1,k2)=V4
      M(k1,k2)=M4
      T(k1,k2)=T4
      p(k1,k2)=p4
      rou(k1,k2)=rou4
      
      k=k2
      ffn(k1)=k2
      ! end successive point on the r-c-line
      
      ! start point D seeking on the r-c-line
      mdotBD=0.0   ! initiate the mass flow rate throught BD on the r-c-line
      muld0=1e6   ! initate muld0 with a huge number
      flagd=0   ! initiate flagd with assume "less than" occurred
      
      if (k1==46) then
         k2=k2+1
         k2=k2-1
      end if
      
      do k2=2,ffn(k1)-1   ! point D start with the second point on the r-c-line,
                          ! end by inverse second point avoids thetad=0
         ! calculate mdotBD	along BD
         arou=(rou(k1,k2-1)+rou(k1,k2))/2
         aV=(V(k1,k2-1)+V(k1,k2))/2
         atheta=(atan(Vv(k1,k2-1)/Vu(k1,k2-1))+atan(Vv(k1,k2)/Vu(k1,k2)))/2
         aalpha=(asin(1/M(k1,k2-1))+asin(1/M(k1,k2)))/2
         ay=(y(k1,k2-1)+y(k1,k2))/2
         mdotBD=mdotBD+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)&
               &*(y(k1,k2-1)-y(k1,k2))
         ! end calculate mdotBD	along BD
         
         ! calculate point along DE
         ! determine the lagrange multiplier at point D
         thetad=atan(Vv(k1,k2)/Vu(k1,k2))
         alphad=asin(1/M(k1,k2))
         C2=V(k1,k2)*cos(thetad-alphad)/cos(alphad)
         C3=y(k1,k2)*rou(k1,k2)*V(k1,k2)**2*(sin(thetad))**2*tan(alphad)
         ! end determine the lagrange multiplier at point D
         
         mdotDE=0.0   ! initiate mass flow rate along DE before the current point
         mdotDE0(1)=0.0   ! initiate mass flow rate along DE before previous point
         k3=k1+1   !   initiate k3 as r-c-line indicator after point D
         mule0=1e6   ! initiate mule0 with a huge number
         
         do   ! point E seeking loop
            y(k3,k2)=y(k3-1,k2)+stepd   !radical location of the point on DE
            ! calculate V and thetae
            call Vtheta(C2,C3,thetad,M(k3-1,k2),y(k3,k2),V(k3,k2),thetae)
            if (abs(thetae)<1e-5) then
               thetae=thetae+1
               thetae=thetae-1
            end if
            ! calculate other properties of the point E
            Vu(k3,k2)=V(k3,k2)*cos(thetae)
            Vv(k3,k2)=V(k3,k2)*sin(thetae)
            call Thermo(V(k3,k2),p(k3,k2),rou(k3,k2),T(k3,k2),atemp,M(k3,k2))
            ! calculate mass flow rate along DE
            arou=(rou(k3-1,k2)+rou(k3,k2))/2
            aV=(V(k3-1,k2)+V(k3,k2))/2
            aVu=(Vu(k3-1,k2)+Vu(k3,k2))/2
            ay=(y(k3-1,k2)+y(k3,k2))/2
            aalpha=(asin(1/M(k3-1,k2))+asin(1/M(k3,k2)))/2
            atheta=(atan(Vv(k3-1,k2)/Vu(k3-1,k2))+atan(Vv(k3,k2)/Vu(k3,k2)))/2
            mdotDE=mdotDE+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha+atheta)*stepd
            ! determine x location of the point E
            x(k3,k2)=x(k3-1,k2)+(y(k3,k2)-y(k3-1,k2))/tan(aalpha+atheta)
            xe=x(k3,k2)   ! xe, Me and pe retaioned for exit the E and D seeking loops
            p0e=p(k3,k2)-0.5*rou(k3,k2)*V(k3,k2)**2*sin(2*thetae)*&
                  &tan(asin(1/M(k3,k2)))
            mule=abs(xe-xed)*abs(p0e-p0)   ! product of error in E loop
            ! if the point is more close to E than previous one
            if (mule < mule0) then
               mule0=mule
            end if
            if (mdotDE > mdotBD) then
               ! interpolate for latest point along DE
               interc=(mdotBD-mdotDE0(k3-k1))/(mdotDE-mdotDE0(k3-k1))
               y(k3,k2)=y(k3-1,k2)+(y(k3,k2)-y(k3-1,k2))*interc
               ! recalculate V and thetae of the latest point
               call Vtheta(C2,C3,thetad,M(k3-1,k2),y(k3,k2),V(k3,k2),thetae)
               ! calculate other properties of the latest point
               Vu(k3,k2)=V(k3,k2)*cos(thetae)
               Vv(k3,k2)=V(k3,k2)*sin(thetae)
               call Thermo(V(k3,k2),p(k3,k2),rou(k3,k2),T(k3,k2),atemp,M(k3,k2))
               ! calculate thrust along DE
               aalpha=(asin(1/M(k3-1,k2))+asin(1/M(k3,k2)))/2
               atheta=(atan(Vv(k3-1,k2)/Vu(k3-1,k2))+atan(Vv(k3,k2)/Vu(k3,k2)))/2
               ! determine x location of the point E
               x(k3,k2)=x(k3-1,k2)+(y(k3,k2)-y(k3-1,k2))/tan(aalpha+atheta)
               xe=x(k3,k2)   ! store the x-value of point E
               p0e=p(k3,k2)-0.5*rou(k3,k2)*V(k3,k2)**2*sin(2*thetae)*&
                     &tan(asin(1/M(k3,k2)))   ! calculate ambient pressure
               ! judge if the design properties obtained
               mule=abs(xe-xed)*abs(p0e-p0)   ! product of error in E loop
               muld=mule   ! store it as muld
               mdotDE0(k3-k1+1)=mdotBD
               exit   ! mass flow rate or xe excceed, exit E loop
            end if
            
            ! if xe > xed and mdotDE <= mdotBD, exit E loop
            if (xe > xed) exit
            ! store mdot before the current point for next interpolation
            mdotDE0(k3-k1+1)=mdotDE
            ! reset thetad for Vtheta subroutine, speed up loops, avoide M(i+1)<M(i)
            thetad=thetae
            k3=k3+1   ! reset k3 for nest E seeking loop
         end do   ! end point E seeking loop
         
         if (muld < muld0) then   !if the error product < previous one
            muld0=muld
            mult=muld
            flagd=0
         else
            flagd=flagd+1
         end if
         
         ! exit loop D when xe > xed
         if (xe > xed) exit
         
      end do
      ! end point D seeking on the r-c-line
      
      ! exit theta loop when xe > xed and p0e < p0
      if ((xe > xed) .and. (p0e < p0)) then   ! all of them crossed, get k4, kd
         if (mult < mult0) then   ! if current error product < previous one, (theta)
            k4=k1   ! this r-c-line is best one
            if (flagd == 0) then   ! last point is best one
               kd=k2
               exit   ! point D obtained, exit
            else   ! previous flagd-th point is best one
               kd=k2-flagd
               exit   ! point D obtained, exit
            end if
         else   ! if current error product >= previous one, (theta)
            k4=k1-1   ! previous r-c-line is best one
            if (flagd0 == 0) then   ! last point is best one
               kd=k20
               exit   ! point D obtained, exit
            else   ! previous flagd-th point is best one
               kd=k20-flagd0
               exit   ! point D obtained, exit
            end if
         end if
      else   ! not all of them crossed
         mult0=mult
      end if
      
      k20=k2   ! revise k20 for next loop
      flagd0=flagd   ! revise flagd0 for next loop
      k1=k1+1	! start thetaB loop loop without point E obtained
   end do   !end thetab loop, inverse wall, expansion wall
   
   ! determine mdotBD, begin
   mdotBD=0.0   ! reset mdotBD
   do i=2,kd
      arou=(rou(k4,i-1)+rou(k4,i))/2
      aV=(V(k4,i-1)+V(k4,i))/2
      ay=(y(k4,i-1)+y(k4,i))/2
      aalpha=(asin(1/M(k4,i-1))+asin(1/M(k4,i)))/2
      atheta=(atan(Vv(k4,i-1)/Vu(k4,i-1))+atan(Vv(k4,i)/Vu(k4,i)))/2
      mdotBD=mdotBD+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)&
            &*(y(k4,i-1)-y(k4,i))
   end do
   ! determine mdotBD, end

   ! DE calculation, begin
   ! determine the lagrange multiplier at point D
   thetad=atan(Vv(k4,kd)/Vu(k4,kd))
   alphad=asin(1/M(k4,kd))
   C2=V(k4,kd)*cos(thetad-alphad)/cos(alphad)
   C3=y(k4,kd)*rou(k4,kd)*V(k4,kd)**2*(sin(thetad))**2*tan(alphad)
   ! end determine the lagrange multiplier at point D
   mdotDE=0.0   ! initiate mass flow rate along DE before the current point
   mdotDE0(1)=0.0   ! initiate mass flow rate along DE before previous point
   k3=k4+1   !   initiate k3 as r-c-line indicator after point D
   k2=kd   ! specify k2 value on the r-c-line
   do
      y(k3,k2)=y(k3-1,k2)+stepd   !radical location of the point on DE
      ! calculate V and thetae
      call Vtheta(C2,C3,thetad,M(k3-1,k2),y(k3,k2),V(k3,k2),thetae)
      if (abs(thetae)<1e-5) then
         thetae=thetae+1
         thetae=thetae-1
      end if
      ! calculate other properties of the point E
      Vu(k3,k2)=V(k3,k2)*cos(thetae)
      Vv(k3,k2)=V(k3,k2)*sin(thetae)
      call Thermo(V(k3,k2),p(k3,k2),rou(k3,k2),T(k3,k2),atemp,M(k3,k2))
      ! calculate mass flow rate along DE
      arou=(rou(k3-1,k2)+rou(k3,k2))/2
      aV=(V(k3-1,k2)+V(k3,k2))/2
      aVu=(Vu(k3-1,k2)+Vu(k3,k2))/2
      ay=(y(k3-1,k2)+y(k3,k2))/2
      aalpha=(asin(1/M(k3-1,k2))+asin(1/M(k3,k2)))/2
      atheta=(atan(Vv(k3-1,k2)/Vu(k3-1,k2))+atan(Vv(k3,k2)/Vu(k3,k2)))/2
      mdotDE=mdotDE+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha+atheta)*stepd
      ! determine x location of the point E
      x(k3,k2)=x(k3-1,k2)+(y(k3,k2)-y(k3-1,k2))/tan(aalpha+atheta)
      xe=x(k3,k2)   ! xe, Me and pe retaioned for exit the E and D seeking loops
      p0e=p(k3,k2)-0.5*rou(k3,k2)*V(k3,k2)**2*sin(2*thetae)*&
            &tan(asin(1/M(k3,k2)))
      ! judge if the design properties obtained
      mule=abs(xe-xed)*abs(p0e-p0)   ! product of error in E loop
      ! mass flow rate interpolation
      if (mdotDE > mdotBD) then
         ! interpolate for latest point along DE
         interc=(mdotBD-mdotDE0(k3-k4))/(mdotDE-mdotDE0(k3-k4))
         y(k3,k2)=y(k3-1,k2)+(y(k3,k2)-y(k3-1,k2))*interc
         ! recalculate V and thetae of the latest point
         call Vtheta(C2,C3,thetad,M(k3-1,k2),y(k3,k2),V(k3,k2),thetae)
         ! calculate other properties of the latest point
         Vu(k3,k2)=V(k3,k2)*cos(thetae)
         Vv(k3,k2)=V(k3,k2)*sin(thetae)
         call Thermo(V(k3,k2),p(k3,k2),rou(k3,k2),T(k3,k2),atemp,M(k3,k2))
         ! calculate thrust along DE
         arou=(rou(k3-1,k2)+rou(k3,k2))/2
         aV=(V(k3-1,k2)+V(k3,k2))/2
         aVu=(Vu(k3-1,k2)+Vu(k3,k2))/2
         ay=(y(k3-1,k2)+y(k3,k2))/2
         aalpha=(asin(1/M(k3-1,k2))+asin(1/M(k3,k2)))/2
         atheta=(atan(Vv(k3-1,k2)/Vu(k3-1,k2))+atan(Vv(k3,k2)/Vu(k3,k2)))/2
         ! determine x location of the point E
         x(k3,k2)=x(k3-1,k2)+(y(k3,k2)-y(k3-1,k2))/tan(aalpha+atheta)
         xe=x(k3,k2)   ! store the x-value of point E
         ye=y(k3,k2)   ! store the y-value of point E
         Me=M(k3,k2)   ! store the Mach number of point E
         p0e=p(k3,k2)-0.5*rou(k3,k2)*V(k3,k2)**2*sin(2*thetae)*&
            &tan(asin(1/M(k3,k2)))   ! calculate ambient pressure
         k5=k3   ! store line number on which E obtained
         mdotDE0(k3-k1+1)=mdotBD
         exit   !exit point E seeking loop with point E obtained
      end if
      
      ! store mdot before the current point for next interpolation
      mdotDE0(k3-k4+1)=mdotDE
      ! reset thetad for Vtheta subroutine, speed up loops, avoide M(i+1)<M(i)
      thetad=thetae
      k3=k3+1   ! reset k3 for nest E seeking loop
   end do
   ! DE calculation, end 
   
!************************************** Region R ***************
   ! calculate mdot along DF, begin
   mdotDF=0.0
   do i=kd+1,ffn(k4)
      arou=(rou(k4,i-1)+rou(k4,i))/2
      aV=(V(k4,i-1)+V(k4,i))/2
      aVu=(Vu(k4,i-1)+Vu(k4,i))/2
      ay=(y(k4,i-1)+y(k4,i))/2
      aalpha=(asin(1/M(k4,i-1))+asin(1/M(k4,i)))/2
      atheta=(atan(Vv(k4,i-1)/Vu(k4,i-1))+atan(Vv(k4,i)/Vu(k4,i)))/2
      mdotDF=mdotDF+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)&
               &*(y(k4,i-1)-y(k4,i))
   end do
   ! calculate mdot along DF, end

   k3=0   ! initiate k3
   do k1=k4+1,k5-1   ! begin c-line loop
      dmdot=mdotDE0(k1-k4+1)+mdotDF   ! get the initial mass flow rate
      do k2=kd-1,1,-1   ! begin point loop on the r-c-line
         if (k2==kinterp) then
            k3=k3+1
            kinterp=0
         end if
         x1=x(k1,k2+1)
         y1=y(k1,k2+1)
         Vu1=Vu(k1,k2+1)
         Vv1=Vv(k1,k2+1)
         x2=x(k1-1,k2-k3)
         y2=y(k1-1,k2-k3)
         Vu2=Vu(k1-1,k2-k3)
         Vv2=Vv(k1-1,k2-k3)
         
         call InnerP()
         if (flagy40==1) then
            flagy40=0
            exit
         end if
         
         ! avoid two y-values close to each other much more, induce dmdot crawling
         if((k2 <= kd-2) .and. (y(k1,k2+1) > y(k1,k2+2)) .and. (&
            &abs(y4-y(k1,k2+1)) < 1e-4)) then   ! interpolation
            k3=k3-1
            kinterp=k2
            call Interp(x(k1,k2+2),y(k1,k2+2),x(k1,k2+1),y(k1,k2+1),x4,y4&
                  &,x(k1,k2),y(k1,k2))
            Vu(k1,k2)=Vu(k1,k2+1)+((x(k1,k2)-x(k1,k2+1))/(x4-x(k1,k2+1)))&
                  &*(Vu4-Vu(k1,k2+1))
            Vv(k1,k2)=Vv(k1,k2+1)+((x(k1,k2)-x(k1,k2+1))/(x4-x(k1,k2+1)))&
                  &*(Vv4-Vv(k1,k2+1))
            V(k1,k2)=sqrt(Vu(k1,k2)**2+Vv(k1,k2)**2)
            call Thermo(V(k1,k2),p(k1,k2),rou(k1,k2),T(k1,k2),atemp,M(k1,k2))
         else   ! normal
            x(k1,k2)=x4
            y(k1,k2)=y4
            Vu(k1,k2)=Vu4
            Vv(k1,k2)=Vv4
            V(k1,k2)=V4
            M(k1,k2)=M4
            T(k1,k2)=T4
            p(k1,k2)=p4
            rou(k1,k2)=rou4
         end if
         
         arou=(rou(k1,k2)+rou(k1,k2+1))/2
         aV=(V(k1,k2)+V(k1,k2+1))/2
         atheta=(atan(Vv(k1,k2)/Vu(k1,k2))+atan(Vv(k1,k2+1)/Vu(k1,k2+1)))/2
         aalpha=(asin(1/M(k1,k2))+asin(1/M(k1,k2+1)))/2
         ay=(y(k1,k2)+y(k1,k2+1))/2
!         mass flow rate before current point
         dmdot0=dmdot
         dmdot=dmdot+2*pi*ay*arou*aV*sin(aalpha)/sin(aalpha-atheta)&
               &*(y(k1,k2)-y(k1,k2+1))
         if (dmdot>=mdot) then   ! begin mass balance
            ffn(k1)=kd-k2+1
            interc=(mdot-dmdot0)/(dmdot-dmdot0)
            x(k1,k2)=x(k1,k2+1)+(x(k1,k2)-x(k1,k2+1))*interc
            y(k1,k2)=y(k1,k2+1)+(y(k1,k2)-y(k1,k2+1))*interc
            Vu(k1,k2)=Vu(k1,k2+1)+(Vu(k1,k2)-Vu(k1,k2+1))*interc
            Vv(k1,k2)=Vv(k1,k2+1)+(Vv(k1,k2)-Vv(k1,k2+1))*interc
            V(k1,k2)=sqrt(Vu(k1,k2)**2+Vv(k1,k2)**2)
            call Thermo(V(k1,k2),p(k1,k2),rou(k1,k2),T(k1,k2),atemp,M(k1,k2))
            exit   ! exit point loop with ending the r-c-line
         end if   ! end mass balance
      end do   ! end point loop on the r-c-line
      k3=k2-1   ! record the translated number for next c-line calculation
      ! data translating
      do i=1,ffn(k1)   
         x(k1,i)=x(k1,i+k3)
         y(k1,i)=y(k1,i+k3)
         Vu(k1,i)=Vu(k1,i+k3)
         Vv(k1,i)=Vv(k1,i+k3)
         V(k1,i)=V(k1,i+k3)
         M(k1,i)=M(k1,i+k3)
         T(k1,i)=T(k1,i+k3)
         p(k1,i)=p(k1,i+k3)
         rou(k1,i)=rou(k1,i+k3)
      end do
      ! end data translating
   end do   ! end c-line loop
!   xe storing
   ffn(k5)=1
   x(k5,1)=x(k5,kd)
   y(k5,1)=y(k5,kd)
   Vu(k5,1)=Vu(k5,kd)
   Vv(k5,1)=Vv(k5,kd)
   V(k5,1)=V(k5,kd)
   M(k5,1)=M(k5,kd)
   T(k5,1)=T(k5,kd)
   p(k5,1)=p(k5,kd)
   rou(k5,1)=rou(k5,kd)

!******************************** performance calculation *********************
   ! Fexit=Fexit act on transonic line + force act on the wall, begin
   Fexit=0.0   ! re-initiate Fexit
   do i=2,Ny   ! thrust act on the transonic line, begin
      arou=(rou(i-1,1)+rou(i,1))/2
      aV=(V(i-1,1)+V(i,1))/2
      aVu=(Vu(i-1,1)+Vu(i,1))/2
      ay=(y(i-1,1)+y(i,1))/2
      aalpha=(asin(1/M(i-1,1))+asin(1/M(i,1)))/2
      atheta=(atan(Vv(i-1,1)/Vu(i-1,1))+atan(Vv(i,1)/Vu(i,1)))/2
      ap=(p(i-1,1)+p(i,1))/2
      Fexit=Fexit+2*pi*ay*(ap-p0e+arou*aV*sin(aalpha)*aVu/sin(aalpha-atheta))&
            &*(y(i,1)-y(i-1,1))
   end do   ! thrust act on the transonic line, end
   
   do i=Ny+1,k5
      call Thrust(p(i,1),y(i,1),M(i,1),p(i-1,1),y(i-1,1),p0e)
   end do
   ! Fexit=Fexit act on transonic line + force act on the wall, end
   
!   Ct=Fexit/(pi*ptemp*yt**2)
!   Ctm=g*((g+1)/2)**(-g/(g-1))*sqrt((g+1)/(g-1))

!***************************** write in file ***************************************
! write the data in the file, begin
   open(unit=110,file='temp.txt')
   ! write in header
   write(110,*) 'MAXT_IRRO Data'
   write(110,*)
   ! write in initial region
   write(110,*) 'Conclusion of the concept'
   write(110,*) 'k4 kd Ae/At L/yt Me pe pa thetae'
   write(110,*) k4,kd,ye**2/yt**2,xe/yt,Me,p(k5,1),p0e,thetae*180/pi
   write(110,*)
   ! write in performance
   write(110,*) 'Nozzle Performance'
   write(110,*) 'mdot mdot1d CDm Fexit F1d Isp Isp1d CDm etaF etaI'
   write(110,*) mdot,mdot1d,CDm,Fexit,F1d,Isp,Isp1d,CDm,etaF,etaI   !,Ct,Ctm
   write(110,*)
   write(110,*) 'Initial Region'
   write(110,*) 'i j x y Vu Vv V M T p rou'
   write(110,*)
   do i=1,Ny
      do j=1,2*i-1
write(110,*) i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),V(i,j),M(i,j),T(i,j),p(i,j),rou(i,j)
      end do
   end do
   ! write in expansion region
   write(110,*)
   write(110,*) 'Expansion wall Region'
   write(110,*) 'i j x y Vu Vv V M T p rou'
   write(110,*)
   do i=Ny+1,k4
      do j=1,ffn(i)
write(110,*) i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),V(i,j),M(i,j),T(i,j),p(i,j),rou(i,j)
      end do
   end do
   ! write in offset region
   write(110,*)
   write(110,*) 'Region R Points'
   write(110,*) 'i j x y Vu Vv V M T p rou'
   write(110,*)
   do i=k4+1,k5
      do j=1,ffn(i)
write(110,*) i,j,x(i,j),y(i,j),Vu(i,j),Vv(i,j),V(i,j),M(i,j),T(i,j),p(i,j),rou(i,j)
      end do
   end do
   close(110)
! write the data in the file, end
   write(*,*) 'complete !'

end program MAXT_IRRO
